      real*8 function amohrc(x,u)
      implicit none
      double precision x(4),u(7),devias(6)
      double precision sigma, sj2, sj3, ermc, fb, tgo, c, cosf,
     1                 sinf, q, sbar, c3ta, theta
      double precision zero, one, two, three, four, six, toler, pi
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1 SIX=6.D0,TOLER=1.0d-10,PI=3.1415926535897D0)

      sigma=zero
      sj2=zero
      sj3=zero
      ermc=zero
      fb=zero
      call getdevias(u,devias,sigma,sj2,sj3)
      tgo=x(1)
      c=x(2)
      cosf=1./dsqrt(1.+tgo*tgo)
      sinf=tgo*cosf
      q=dsqrt(three*sj2)
      sbar=dsqrt(sj2)
      if(dabs(sbar).lt.toler) then
        fb=-1.0d0*c
        amohrc=fb
        return
      endif
      c3ta=1.5d0*dsqrt(three)*sj3/(sbar**3)
      if(c3ta.gt.1.0d0) then
      	c3ta=1.0d0
      else if(c3ta.lt.-1.0d0) then
      	c3ta=-1.0d0
      endif
      theta=dacos(c3ta)/three
      if(theta.le.toler) then
      	theta=zero
      elseif(theta.gt.pi/three) then
      	theta=pi/three
      endif

      call getrmc(theta,tgo,cosf,ermc)
      fb=ermc*q+sigma*tgo-c
      amohrc=fb
      return
      end

      subroutine getrmc(theta,tgo,cosf,ermc)
      implicit none
      double precision theta, tgo, cosf, ermc
      double precision sinf, sina, cosa, sinb, cosb
      double precision sinab, cosab
      double precision zero, one, two, three, four,
     1                 six, toler, pi
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1 SIX=6.D0,TOLER=1.0d-10,PI=3.1415926535897D0)
c      print *,theta,tgo,cosf,'--theta,tgo,cosf--'
      sinf=dsqrt(1.0d0-cosf*cosf)
      sina=dsin(theta)
      cosa=dcos(theta)
      sinb=dsin(pi/three)
      cosb=dcos(pi/three)
      sinb=dsqrt(three)/two
      cosb=one/two
      sinab=sina*cosb+cosa*sinb
      cosab=cosa*cosb-sina*sinb
      ermc=sinab/(dsqrt(three)*cosf)+cosab*tgo/three
c      print *,sinab,cosf,cosab,tgo,'---sinab,cof,csab,tgo,ermc--'
      return
      end

      subroutine getrmw(theta,tgo,cosf,ermw)
      implicit none
      double precision theta, tgo, cosf, ermw
      double precision ermc, theta1, cosa, sinf, e, s, ss0, ss1, ss
      double precision zero, one, two, three, four, five, toler, pi
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1          FIVE=5.d0, TOLER=1.0d-10,PI=3.1415926535897D0)
      ermc=zero
      theta1=pi/three
      call getrmc(theta1,tgo,cosf,ermc)
      sinf=tgo*cosf
      ermc=three-sinf
      ermc=ermc/cosf/6.0d0
      cosa=dcos(theta)
c      sinf=dsqrt(1.0d0-cosf*cosf)
      e=(three-sinf)/(three+sinf)
      e=0.6d0
      s=four*(one-e**2)*cosa*cosa+(two*e-one)**2
      ss0=two*(one-e*e)*cosa
      ss1=dsqrt(four*(one-e*e)*cosa*cosa+five*e*e-four*e)
      ss1=(two*e-one)*ss1
      ss=ss0+ss1
      ermw=s*ermc/ss
      return
      end

      subroutine gtheta(u,theta)
      implicit none
      double precision u(7),devias(6), theta
      double precision sigma, sj2, sj3, sbar, q, c3ta
      double precision zero, one, two, three, toler, pi
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,
     1          TOLER=1.0d-10,PI=3.1415926535897D0)
      sigma=zero
      sj2=zero
      sj3=zero
      theta=zero
      call getdevias(u,devias,sigma,sj2,sj3)
      sbar=dsqrt(sj2)
      q=dsqrt(three*sj2)
      if(abs(sbar).lt.toler) then
      c3ta=1.0d0
      theta=zero
      return
      endif
      c3ta=1.5d0*dsqrt(three)*sj3/(sbar**3)
      if(c3ta.gt.1.0d0) then
      	c3ta=1.0d0
      else if(c3ta.lt.-1.0d0) then
      	c3ta=-1.0d0
      endif
      theta=dacos(c3ta)/three
      if(theta.le.toler) then
      	theta=zero
      elseif(theta.gt.pi/three) then
      	theta=pi/three
      endif
      return
      end

      subroutine getg(x,u,fg)
      implicit none
      double precision x(4),u(7),devias(6)
      double precision sigma, sj2, sj3, theta, ermw
      double precision q, tgo, c, cosf, sinf, pusi, tgp
      double precision ec1, ec2, eg, fg
      double precision zero, one, two, three, four,
     !                 six, toler, pi, em
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1 SIX=6.D0,TOLER=1.0d-10,PI=3.1415926535897D0,em=0.1d0)
      sigma=zero
      sj2=zero
      sj3=zero
      theta=zero
      ermw=zero
      call getdevias(u,devias,sigma,sj2,sj3)
      q=dsqrt(three*sj2)
      tgo=x(1)
      c=x(2)
      cosf=1./dsqrt(1.+tgo*tgo)
      sinf = tgo*cosf
C      sinf=dsqrt(1.0d0-cosf*cosf)
      pusi=x(4)
      tgp=tand(pusi)
      call gtheta(u,theta)
      call getrmw(theta,tgo,cosf,ermw)
      sigma=(u(1)+u(2)+u(3))/three
      ec1=em*c*tgp
      ec2=ermw*q
      eg=ec1*ec1+ec2*ec2
      fg=dsqrt(eg)+sigma*tgp
      return
      end

      subroutine  getdevias(u,s,d,sj2,sj3)
      implicit none
      double precision u(7),s(6), d, sj2, sj3
      double precision uk, sj1
      d=(u(1)+u(2)+u(3))/3.0d0
      s(1)=u(1)-d
      s(2)=u(2)-d
      s(3)=u(3)-d
      s(4)=u(4)
      s(5)=u(5)
      s(6)=u(6)
      uk=u(7)
      sj1=s(1)+s(2)+s(3)
      sj2=(s(1)**2+s(2)**2+s(3)**2)/2.0d0+s(4)**2
      sj2=sj2+s(5)**2+s(6)**2
      sj3=s(1)*s(2)*s(3)+2.0d0*s(4)*s(5)*s(6)-s(1)*s(4)*s(4)
      sj3=sj3-s(2)*s(5)*s(5)-s(3)*s(6)*s(6)
      return
      end

      subroutine getyield(x,u,d,shear,ha12,ddf,kq,yield,ialpha)
      implicit none
      integer kq, ialpha
      double precision shear, ha12
      double precision x(4),u(7),v(7),h(7),df(7),ddf(6),d(3,3)
      double precision dg(7),ddg(6)
      double precision theta, fg2, fg1, theta1, dtheta, f2, f1
      double precision c, a, am, h11, h12, h22, cg
      integer i, j, k
      double precision zero, three, toler, pi
      PARAMETER(ZERO=0.D0,THREE=3.D0,TOLER=1.0d-10,PI=3.1415926535897D0)
      double precision yield
      external yield

      theta=zero
      fg2=zero
      fg1=zero
      do k=1,7
       v(k)=u(k)
      end do
      call gtheta(v,theta)
      theta1=pi/three
      dtheta=theta-theta1
      if(theta.le.toler) then
      	call getdf0(x,v,df,dg)
      	goto 1
      elseif(dabs(dtheta).lt.toler) then
      	call getdf1(x,v,df,dg)
        goto 1
      endif
      do k=1,6
       h(k) = 1.0d-3
      end do
      if(kq.eq.1) h(7)=1.0e-3
      if(kq.eq.2) h(7)=1.0e-6
      if(kq.eq.3) h(7)=1.0e-6
      do k=1,7
       v(k)=u(k)+h(k)
       f2=yield(x,v)
       call getg(x,v,fg2)
       v(k)=u(k)-h(k)
       call getg(x,v,fg1)
       f1=yield(x,v)
       df(k)=(f2-f1)/h(k)*0.5d0
       dg(k)=(fg2-fg1)/h(k)*0.5d0
       v(k)=u(k)
      end do
1     continue
      do i=1,7
      	dg(i)=df(i)
      enddo
      do i=4,6
       ddf(i)=shear*df(i)
       ddg(i)=shear*dg(i)
      end do
      do i=1,3
       c=0.0d0
       cg=0.0d0
       do j=1,3
        c=c+d(i,j)*df(j)
        cg=cg+d(i,j)*dg(j)
       end do
       ddf(i)=c
       ddg(i)=cg
      end do
      a=0.0d0
      am=0.0d0
      h11=0.0d0
      h12=0.0d0
      h22=0.0d0
      ha12=0.0d0
      do i=1,6
       a=a+df(i)*ddf(i)
       h11=h11+df(i)*ddf(i)
       h12=h12+df(i)*ddg(i)
       h22=h22+dg(i)*ddg(i)
       if (kq.eq.1) am=am+u(i)*dg(i)
       if ((kq.eq.2).and.(i.le.3)) am=am+dg(i)
       if (kq.eq.3) am=am+dg(i)*dg(i)
      end do
      if (kq.eq.3) am=sqrt(am)
      am=-am*df(7)
      if(ialpha.eq.0)then
        ha12=h12+am
        a=a+am
      else
c       if (am.ge.-a) then
c        a=a+am+a*3.5
c       else
c        a=a+am-a*2.0
c       end if
c      end if
       if (am.ge.-h12) then
        ha12=h12+am+h12*3.5d0
       else
        ha12=h12+am-h12*2.0d0
       end if
      end if
C      print *,'In afg: ddf =', ddf
      return
      end

      subroutine getdf0(x,v,df,dg)
      implicit none
      double precision x(4),v(7),df(7),dg(7),devias(6)
      double precision sigma, sj2, sj3, sbar, tgo, c, pusi,
     1                 cosf, sinf, e, c1, ermc, c2, cg1,
     2                 ermw, cg2
      integer          i
      double precision zero, one, two, three, four, six,
     1                 toler, pi
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1 SIX=6.D0,TOLER=1.0d-10,PI=3.14159276535897D0)
      sigma=zero
      sj2=zero
      sj3=zero
      call getdevias(v,devias,sigma,sj2,sj3)
      sbar=dsqrt(sj2)
      if(sbar.lt.toler) then
      do i=1,7
      	df(i)=0.0d0
      	dg(i)=0.0d0
      enddo
      return
      endif
      tgo=x(1)
      c=x(2)
      pusi=x(4)
      cosf=1./dsqrt(1.+tgo*tgo)
      sinf=dsqrt(1.0d0-cosf*cosf)
      e=(three-sinf)/(three+sinf)
      c1=tgo
      ermc=1.0d0/(2.0d0*cosf)+tgo/6.0d0
      c2=dsqrt(three)*ermc
      df(1)=C1/THREE+C2*DEVIAS(1)/(two*SBAR)
      df(2)=C1/THREE+C2*DEVIAS(2)/(two*SBAR)
      df(3)=C1/THREE+C2*DEVIAS(3)/(two*SBAR)
      df(4)=c2*devias(4)/sbar
      df(5)=c2*devias(5)/sbar
      df(6)=c2*devias(6)/sbar
      df(7)=zero
      cg1=tand(pusi)
      ermw=1.0/(2.0d0*cosf)-tgo/6.0d0
      cg2=dsqrt(three)*e*ermw
      dg(1)=Cg1/THREE+Cg2*DEVIAS(1)/(two*SBAR)
      dg(2)=Cg1/THREE+Cg2*DEVIAS(2)/(two*SBAR)
      dg(3)=Cg1/THREE+Cg2*DEVIAS(3)/(two*SBAR)
      dg(4)=cg2*devias(4)/sbar
      dg(5)=cg2*devias(5)/sbar
      dg(6)=cg2*devias(6)/sbar
      dg(7)=zero
      return
      end

      subroutine getdf1(x,v,df,dg)
      implicit none
      double precision x(4),v(7),df(7),dg(7),devias(6)
      double precision sigma, sj2, sj3, sbar, tgo, c, pusi, cosf,
     1                 sinf, e, c1, ermc, c2, cg1, ermw, cg2
      integer          i
      double precision zero, one, two, three, four, six,
     1                 toler, pi
      PARAMETER(ZERO=0.D0,ONE=1.D0,TWO=2.D0,THREE=3.D0,FOUR=4.D0,
     1 SIX=6.D0,TOLER=1.0d-10,PI=3.1415922535897D0)
      sigma=zero
      sj2=zero
      sj3=zero
      call getdevias(v,devias,sigma,sj2,sj3)
      sbar=dsqrt(sj2)
      if(sbar.lt.toler) then
      do i=1,7
      	df(i)=0.0d0
      	dg(i)=0.0d0
      enddo
      return
      endif
      tgo=x(1)
      c=x(2)
      pusi=x(4)
      cosf=1.d0/dsqrt(1.+tgo*tgo)
      sinf=tgo*cosf
      e=(three-sinf)/(three+sinf)
      c1=tgo
      ermc=1.0/(2.0d0*cosf)-tgo/6.0d0
      c2=dsqrt(three)*ermc
      df(1)=C1/THREE+C2*DEVIAS(1)/(two*SBAR)
      df(2)=C1/THREE+C2*DEVIAS(2)/(two*SBAR)
      df(3)=C1/THREE+C2*DEVIAS(3)/(two*SBAR)
      df(4)=c2*devias(4)/sbar
      df(5)=c2*devias(5)/sbar
      df(6)=c2*devias(6)/sbar
      df(7)=zero
      cg1=tand(pusi)
      ermw=1.0/(2.0d0*cosf)-tgo/6.0d0
      cg2=dsqrt(three)*ermw
      dg(1)=Cg1/THREE+Cg2*DEVIAS(1)/(two*SBAR)
      dg(2)=Cg1/THREE+Cg2*DEVIAS(2)/(two*SBAR)
      dg(3)=Cg1/THREE+Cg2*DEVIAS(3)/(two*SBAR)
      dg(4)=cg2*devias(4)/sbar
      dg(5)=cg2*devias(5)/sbar
      dg(6)=cg2*devias(6)/sbar
      dg(7)=zero
      return
      end